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Abstract. We propose here a number of approaches to implement con- 
straint propagation for arithmetic constraints on integer intervals. To 
this end we introduce integer interval arithmetic. Each approach is ex- 
plained using appropriate proof rules that reduce the variable domains. 
We compare these approaches using a set of benchmarks. 



1 Preliminaries 
1.1 Introduction 

The subject of arithmetic constraints on reals has attracted a great deal of atten- 
tion in the literature. For some reason arithmetic constraints on integer intervals 
have not been studied even though they are supported in a number of constraint 
programming systems. In fact, constraint propagation for them is present in 
ECL*PS% SICStus Prolog, GNU Prolog, ILOG Solver and undoubtedly most of 
the systems that support constraint propagation for linear constraints on integer 
intervals. Yet, in contrast to the case of linear constraints — see notably [5] — 
we did not encounter in the literature any analysis of this form of constraint 
propagation. 

In this paper we study these constraints in a systematic way. It turns out 
that in contrast to linear constraints on integer intervals there are a number of 
natural approaches to constraint propagation for these constraints. 

To define them we introduce integer interval arithmetic that is modeled after 
the real interval arithmetic, see e.g., [6]. There are, however, essential differences 
since we deal with integers instead of reals. For example, multiplication of two 
integer intervals does not need to be an integer interval. In passing by we show 
that using integer interval arithmetic we can also define succinctly the well- 
known constraint propagation for linear constraints on integer intervals. In the 
second part of the paper we compare the proposed approaches by means of a set 
of benchmarks. 



1.2 Constraint Satisfaction Problems 

We review here the standard concepts of a constraint and of a constraint satis- 
faction problem. Consider a sequence of variables X := xi, . . .,Xn where n > 0, 



with respective domains Di, . . ., D„ associated with them. So each variable Xi 
ranges over the domain Di. By a constraint C on X we mean a subset of 
Di X . . . X Dn- Given an element d := di, . . .,dn of Di x . . . x D„ and a sub- 
sequence Y := Xi^,. . ., Xii, of X we denote by d\Y] the sequence , . . ., di^ . In 
particular, for a variable Xi from X, d[xi\ denotes di. 

A constraint satisfaction problem, in short CSP, consists of a finite se- 
quence of variables X with respective domains V, together with a finite set C of 
constraints, each on a subsequence of X. We write it as (C ; x\ & Di, . . .,Xn & 
Z)„), where X :~ xi, . . ., .t„ and T) := Di, . . ., _D„. 

By a solution to (C ; xi G Di, . . ., a;„ G Z)„) wc mean an clement d, G -Di x 
... X Dn such that for each constraint C G C on a sequence of variables X we have 
d[X] G C. We call a CSP consistent if it has a solution and inconsistent if it 
does not. Two CSPs with the same sequence of variables are called equivalent 
if they have the same set of solutions. In what follows we consider CSPs the 
constraints of which are defined in a simple language and identify the syntactic 
description of a constraint with its meaning being the set of tuples that satisfy 
it. 

We view constraint propagation as a process of transforming CSPs that 
maintains their equivalence. In what follows we define this process by means of 
proof rules that act of CSPs and preserve equivalence. An interested reader can 
consult [1] for a precise explanation of this approach to describing constraint 
propagation. 

1.3 Arithmetic Constraints 

To define the arithmetic constraints use the alphabet that comprises 

— variables, 

— two constants, and 1, 

— the unary minus function symbol '— ', 

— three binary function symbols, '+','— 'and all written in the infix notation. 

By an arithmetic expression we mean a term formed in this alphabet and 
by an arithmetic constraint a formula of the form 

s opt, 

where s and t are arithmetic expressions and op G {<, <, =, 7^, >, >}. For exam- 
ple 

x^ -y"^ ■ + 3x-y^ ■ <10 + 4x^ ■ ■ 7? - ■ x^ ■ z'' (1) 

is an arithmetic constraint. Here x^ is an abbreviation for and 
similarly with the other expressions. If '•' is not used in an arithmetic constraint, 
wc call it a linear constraint. 

By an extended arithmetic expression we mean a term formed in the 
above alphabet extended by the unary function symbols and ' for each 



n > 1 and the binary function symbol '/' written in the infix notation. For 
example 

. z4)/(,^2 . „5) (2) 

is an extended arithmetic expression. Here, in contrast to the above is a term 
obtained by applying the function symbol to the variable x. The extended 
arithmetic expressions will be used only to define constraint propagation for the 
arithmetic constraints. 

Fix now some arbitrary linear ordering -< on the variables of the language. 
By a monomial we mean an integer or a term of the form 

where A; > 0, different variables ordered w.r.t. and a is a non- 

zero integer and ni, . . .,nk are positive integers. We call then x^'^ • . . . • x^'' the 
power product of this monomial. 

Next, by a polynomial we mean a term of the form 

where n > 0, at most one monomial is an integer, and the power products 
of the monomials mi, . . .,m„ arc pairwisc different. Finally, by a polynomial 
constraint we mean an arithmetic constraint of the form s op b, where s is a 
polynomial with no monomial being an integer, op G {<, <, =, ^, >, >}, and b is 
an integer. It is clear that by means of appropriate transformation rules we can 
transform each arithmetic constraint to a polynomial constraint. For example, 
assuming the ordering a; ^ y ^ ^ on the variables, the arithmetic constraint (1) 
can be transformed to the polynomial constraint 

2x^ ■y'^-z'^- 4x^ -y^ ■ z'^ + 3x-y^ ■ <10 

So, without loss of generality, from now on we shall limit our attention to the 
polynomial constraints. 

Next, let us discuss the domains over which we interpret the arithmetic con- 
straints. By an integer interval, or an interval in short, we mean an expres- 
sion of the form 

[a..b] 

where a and b are integers; [a..b] denotes the set of all integers between a and b, 
including a and fe. If a > b, we call [a..b] the empty interval and denote it by 
0. Finally, by a range we mean an expression of the form 

X e I 

where a; is a variable and I is an interval. 
2 Integer Set Arithmetic 

To reason about the arithmetic constraints we employ a generalization of the 
arithmetic operations to the sets of integers. 



2.1 Definitions 



For X, Y sets of integers we define the following operations: 



— addition: 



X + Y 



{x + y\xeX,yeY} 



— subtraction: 



X-Y 



{x-y\xGX,yGY} 



— multiplication: 



X-Y 



{x-y\x&X,y&Y}, 



— division: 



X/Y ■.= {ueZ\3xe X3y eYu-y = x} 



— exponentiation: 



:= {x"-\xG X} 



for each natural number n > 0, 
— root extraction: 



\/X := {a; e ^ I a;" e X}, 



for each natural number n > 0. 

All the operations except division are defined in the expected way. We shall 
return to it at the end of Section 6. At the moment it suffices to note the division 
operation is defined for all sets of integers, including F = and Y = {0}. This 
division operation corresponds to the following division operation on the sets of 
reals introduced in [8]: 



For a (n integer or real) number a and op € {+, —,•)/} we identify a op X with 
{a} op X and X opa with X op {a}. 

To present the rules we are interested in we shall also use the addition and 
division operations on the sets of real numbers. Addition is defined in the same 
way as for the sets of integers, and division is defined above. In [6] it is explained 
how to implement these operations. 

Further, given a set A of integers or reals, we define 



-A := {x G Z \ 3a G Ax > a}. 

When limiting our attention to intervals of integers the following simple ob- 
servation is of importance. 

Note 1. For X,Y integer intervals and a an integer the following holds: 




< 



A := {x G Z \ 3a G Ax < a} 



— X nY, X + Y,X — Y are integer intervals. 



— X/{a} is an integer interval. 

— X - Y does not have to be an integer interval, even li X = {a} or y = {o}. 

— X/Y does not have to be an integer interval. 

— For each n > 1 X" does not have to be an integer interval. 

— For odd n > 1 \/X is an integer interval. 

— For even n > 1 \/X is an integer interval or a disjoint union of two integer 

intervals. □ 

For example we have 

[2..4] + [3..8] = [5..12], 

[3..7] - [1..8] = [-5..6], 

[3..3]-[1..2] = {3,6}, 
[3..5]/[-1..2] = {-5,-4,-3,2,3,4,5}, 
[-3..5]/[-1..2] = Z, 
[1..2f = {1,4}, 



\/[-30..100] ^ [-3. .4] 



^[-100..9] = [-3..3], 
^/[T::9] = [-3..-l]u[1..3]. 

To deal with the problem that non-interval domains can be produced by some 
of the operations we introduce the following operation on the subsets of the set 
of the integers Z: 



\ smallest integer interval containing X if X is finite, 
\ Z otherwise. 



int{X) : 

For example mf([3..5]/[-1..2]) = [-5. .5] and mt([-3..5]/[-1..2]) = Z. 



2.2 Implementation 

To define constraint propagation for the arithmetic constraints on integer in- 
tervals we shall use the integer set arithmetic, mainly limited to the integer 
intervals. This brings us to the discussion of how to implement the introduced 
operations on the integer intervals. Since we are only interested in maintaining 
the property that the sets remain integer intervals or the set of integers Z we 
shall clarify how to implement the intersection, addition, subtraction and root 
extraction operations of the integer intervals and the int{.) closure of the mul- 
tiplication, division and exponentiation operations on the integer intervals. The 
case when one of the intervals is empty is easy to deal with. So we assume that 
we deal with non-empty intervals [a..b] and [c.rf], that is a <b and c< d. 



Intersection, addition and subtraction It is easy to see that 

[a..b] n [c.d] = [max{a,c)..min{b,d)], 

[a..b] + [c.d] = [a + c .. b + d], 

[a..b] — [c.d] = [a — d .. b — c]. 

So the interval intersection, addition, and subtraction are straightforward to 
implement. 

Root extraction The outcome of the root extraction operator applied to an inte- 
ger interval will be an integer interval or a disjoint union of two integer intervals. 
Wc shall explain in Section 4 why it is advantageous not to apply int{.) to the 
outcome. This operator can be implemented by means of the following case 
analysis. 



Case 1. Suppose n is odd. Then 



Case 2. Suppose n is even and 6 < 0. Then 



rb 



vTm = 0- 

Case 3. Suppose n is even and 6 > 0. Then 

^/Ml = [- [l ^ij •• - [l ^ij] u [[l ^[] .. [l ^ij] 

where a+ := max{Q,a). 



Multiplication For the remaining operations we only need to explain how to 
implement the int{.) closure of the outcome. First note that 

int{[a..b] ■ [c..d\) = [min{A)..max{A)], 

where A= {a ■ c,a ■ d,b ■ c,b ■ d}. 

Using an appropriate case analysis wc can actually compute the bounds of 
int{[a..b] ■ [c.d]) directly in terms of the bounds of the constituent intervals. 

Division In contrast, the int{.) closure of the interval division is not so straight- 
forward to compute. The reason is that, as wc shall sec in a moment, wc cannot 
express the result in terms of some simple operations on the interval bounds. 

Consider non-empty integer intervals [a..b] and [c.d]. In analyzing the out- 
come of int{[a..b]/[c..d]) we distinguish the following cases. 

Case 1. Suppose e [a..b] and € [c.d\. 

Then by definition int{[a..b\/[c..d\) = Z. For example, 



mi([-1..100]/[-2..8]) = Z. 



Case 2. Suppose ^ [a..h] and c = d = 0. 

Then by definition int{[a..h]/[c..d\) = 0. For example, 



mt([10..100]/[0..0]) = 0. 

Case 3. Suppose ^ [a..b] and c < and < d. 
It is easy to see that then 

int{[a..b]/[c..d\) = [— e..e], 

where e = max{\a\, \b\). For example, 

mi([-100.. - 10]/[-2..5]) = [-100. .100]. 

Case 4- Suppose ^ [a..b] and cither c = and c? 7^ or c ^ and d = 0. 
Then int{[a..b]/[c..d]) = int{[a..b]/{[c..d] - {0})). For example 

mi([1..100]/[-7..0]) = mt([1..100]/[-7.. - 1]). 

This allows us to reduce this case to Case 5 below. 
Case 5. Suppose ^ [c.d]. 

This is the only case when we need to compute int{[a..b]/[c..d]) indirectly. 
First, observe that we have 

int{[a..b]/[c..d]) C [\min{A)] .. lmax{A)\], 

where A = {a/c, a/d, b/c, b/d}. 

However, the equality does not need to hold here. Indeed, note for example 
that mt([155..161]/[9..11]) = [16.. 16], whereas for A = {155/9,155/11,161/9, 
161/11} we have \min{A)~\ = 15 and [maa;(j4)J = 17. The problem is that the 
value 16 is obtained by dividing 160 by 10 and none of these two values is an 
interval bound. 

This complication can be solved by preprocessing the interval [c.d] so that 
its bounds are actual divisors of an element of [a..b]. First, we look for the least 
c' e [c.rf] such that 3x G [a..b] 3u G Z u-c' = x. Using a case analysis, the latter 
property can be established without search. Suppose for example that a > and 
c > 0. In this case, if c' ■ [^J > a, then c' has the required property. Similarly, 
we look for the largest d' G [c.d] for which an analogous condition holds. Now 
int{[a..b]/[c.d]) = [\min{A)]..lmax{A)\], where A = {a/c', a/d', b/c', b/d'}. 

Exponentiation The int{.) closure of the interval exponentiation is straightfor- 
ward to implement by distinguishing the following cases. 

Case 1. Suppose n is odd. Then 

int{[a..bf) = [a".. 6"]. 
Case 2. Suppose n is even and < a. Then 



int{[a..bY-) = [a".. 6"]. 



Case 3. Suppose n is even and 6 < 0. Then 

mt([a..6]") = [6".. a"]. 
Case 4- Suppose n is even and a < and < 6. Then 
mi([a..6]") = [0..maa;(a", 6")]. 

2.3 Correctness Lemma 

Given now an extended arithmetic expression s each variable of which ranges 
over an integer interval, we define int{s) as the integer interval or the set Z 
obtained by systematically replacing each function symbol by the application of 
the int{.) operation to the corresponding integer set operation. For example, for 
the extended arithmetic expression s := -^(y^ • z^)/{x'^ ■ u^) of (2) we have 

int(s) = int{^int{int{Y'^) ■ int{Z^))/int(int{X'^) ■ int{U^))), 

where x ranges over X, etc. 

The discussion in the previous subsection shows how to compute int{s) given 
an extended arithmetic expression s and the integer interval domains of its vari- 
ables. 

The following lemma is crucial for our considerations. It is a counterpart 
of the so-called 'Fundamental Theorem of Interval Arithmetic' established in 
[7]. Because we deal here with the integer domains an additional assumption is 
needed to establish the desired conclusion. 

Lemma 1 (Correctness). Let s he an extended arithmetic expression with the 
variables Xi,...,a;„. Assume that each variable Xi of s ranges over an integer 
interval Xi. Choose ai € X^ for i € [l..n] and denote by s{ai, . . .,a„) the result 
of replacing in s each occurrence of a variable xi by ai . 

Suppose that each subexpression ofs{ai, . . .,a„) evaluates to an integer. Then 
the result of evaluating s(ai, . . .,a„) is an element of int{s). 

Proof. The proof follows by a straightforward induction on the structure of s. 

□ 

3 An Intermezzo: Constraint Propagation for Linear 
Constraints 

Even though we focus here on arithmetic constraints on integer intervals, it is 

helpful to realize that the integer interval arithmetic is also useful to define 
in a succinct way the well-known rules for constraint propagation for linear 
constraints. To this end consider first a constraint Sf^-^ai ■ Xi = b, where n > 0, 
ai,...,a„ arc non-zero integers, xi,...,Xn are different variables, and b is an 
integer. To reason about it we can use the following rule parametrized by j € 
[l..n]: 



where 
— for i ^ j 



LINEAR EQUALITY 
{Uj^iai ■Xi = b ; xi € Di,. . .,Xn € Dn) 



D' := Di 



D'j := Dj n int({b - Sie[i..n]-{j}a'i ■ Xi)/a^. 

Note that by virtue of Note 1 

D'j = Dj n (6 - i:j£[i..„]_{j}mi(ai • Di))/aj. 

To see that this rule preserves equivalence suppose that for some di € 
Di, . . .,dn & Dn we have Si^iUi ■ di = b. Then for j e [l..n] we have 

dj = (& - Si(z[i,,n]-{j}ai ■ di)/aj 
which by the Correctness Lemma 1 implies that 

dj e int{{b - i7j£[i..„]_{j}ai • Xi)/aj^ , 

i.e., dj e D'j. 

Next, consider a constraint S"_iai ■ Xi < b, where ai, . . ., o„, Xi, . . .,x„ and b 
are as above. To reason about it we can use the following rule parametrized by 
j e [l..n]: 

LINEAR INEQUALITY 

(rjLiO,, • <b; xiGDu..., ,t„ G D„) 



where 
— for i ^ j 



{^i=iaz ■Xi<b; xi € D[,...,Xn e D'J 



A' := A 



-D^- := -Dj n {-int(b - Z'ig[i..„]_{ji,ai • Xi)/aj) 

To see that this rule preserves equivalence suppose that for some di € 
r>i, . . .,dn S Dn we have i^-Li^Sj ■ di < b. Then ■ dj < b — 'S'i£[i..n]-{j}<ii " 
By the Correctness Lemma 1 

^ ~ ^ie[i..n]-{j}0'i ■ di e mi(6 - -S'i£[i..n]-{j}ai • a;,), 
so by definition 

aj ■ dj e- int{b - Si^[i__n]-{j}ai ■ Xi) 

and consequently 

dj e- int{b- Si^[i..n]-{j}ai ■ Xi)/aj 
This implies that dj G D'j. 



4 Constraint Propagation: First Approach 



We now move on to a discussion of constraint propagation for the arithmetic 
constraints on integer intervals. To illustrate the first approach consider the 
following example. Consider the constraint 

x^y - a; < 40 

and the ranges x G [1..100] and y G [1..100]. We can rewrite it as 



X < 



40- 



(3) 



since x assiimes integer values. The maximum value the expression on the right- 
hand side can take is [v^l40j , so we conclude x < 5. By reusing (3), now with the 
information that x € [1..5], we conclude that the maximum value the expression 
on the right-hand side of (3) can take is actually [v^isj, from which it follows 
that .X < 3. 

In the case of y we can isolate it by rewriting the original constraint as 
from which it follows that y < 41, since by assumption x > 1. So 



y 



< 



40 



we could reduce the domain of x to [1..3] and the domain of y to [1..41]. This 
interval reduction is optimal, since x = 1,2/ = 41 and x = 3,y = 1 are both 
solutions to the original constraint x^y — x < 40. 

More formally, we consider a polynomial constraint S^inii = b where m > 0, 
no monomial rrii is an integer, the power products of the monomials are pairwise 
different and b is an integer. Suppose that xi,. . .,Xn are its variables ordered 
w.r.t. -<. 

Select a non-integer monomial me and assume it is of the form 



Vk ' 



where A; > 0, yi 
integer and m, 
variable in {xi, . 
proof rule: 



where 
— for i ^ j 



,.. .,yk are different variables ordered w.r.t. ^, a is a non-zero 

..,nfc are positive integers. So each yi variable equals to some 
. .,x„}. Suppose that yp equals to xj. We introduce the following 

POLYNOMIAL EQUALITY 

(r-Limj =b ; a;i g £>i, . . .,a;„ g £>„) 
(rrLimi =b ; xiG D[,...,XnG D'„) 



D' := A, 



and 



Dj := int (^D^ D int ((6 - Sie[i..m]-{e}'mi)/s) 



s .— a y^ ... ... yf. . 



To see that this rule preserves equivalence choose some rfi G Di, . . ., rf„ e Z)„. 
To simplify the notation, given an extended arithmetic expression t denote by t' 
the result of evaluating t after each occurrence of a variable Xi is replaced by d,. 

Suppose that S^-^^m'^ = b. Then 

(f/ ■ S = b- i^i6[l..ml-{0™i' 

SO by the Correctness Lemma 1 applied to 6 — i^ie[i..TO]-{«}"^i a^iid to s 
d"" e int{b- Si^^i„rn\-{e}tni)lint{s). 

Hence 

'^ie[l..m]-{£}™i)/^'^*(*) 

and consequently 



dj e int y:)j n y ((6 - i^ie[i..m]- {£}'Tii)/s) 

i.e., e l?^-. 



Note that we do not apply int{.) to the outcome of the root extraction op- 
eration. For even Up this means that the second operand of the intersection can 
be a union of two intervals, instead of a single interval. To see why this is de- 
sirable, consider the constraint — y = in the presence of ranges x G [0..10], 
y e [25.. 100]. Using the int{.) closure of the root extraction we would not be 
able to update the lower bound of x to 5. 

Next, consider a polynomial constraint S'^-^rrii < b. Below we adopt the 
assumptions and notation used when defining the POLYNOMIAL EQUALITY 
rule. To formulate the appropriate rule we stipulate that for extended arithmetic 
expressions s and t 

int{{^s)/t) := ^Qn ^Q, 

with Q = {^int{s))/int{t). 

To reason about this constraint we use the following rule: 

POLYNOMIAL INEQUALITY 

{SjL-imi <b ; xi e Di, . . .,Xn & Dn) 
{S^rrii <b; xi€D[,.. .,Xn e D'„) 

where 
— for i ^ j 

D', := A, 




To prove that this rule preserves equivalence choose some di G Di , . . . , c?„ G 
Dn- As above given an extended arithmetic expression t we denote by t' the 
result of evaluating t when each occurrence of a variable Xi in t is replaced by 
di. 

Suppose that S^iin'i < b. Then 

d^P . s' <b- i^ie[l..m]-{<!}W-- 

By the Correctness Lemma 1 

b - ^ie[i..m]-{«}K e int{b - Sieii..m]-{e}'mi), 

so by definition 

rf"" • s' e- int{b- Si^[i„m]-{e}mi). 
Hence by the definition of the division operation on the sets of integers 

d"" G- int{b- IJi^[i„m]-{e}mi)/int{s) 

Consequently 

dj G ^int{b - rig[i..„]_|^|mi)/ini(s) 

This implies that G Dj. 

Note that the set -int{b— S^^^i ^]_^gjmi) is not an interval. So to properly 
implement this rule we need to extend the implementation of the division oper- 
ation discussed in Subsection 2.2 to the case when the numerator is an extended 
interval. Our implementation takes care of this. 

In an optimized version of this approach we simplify the fractions of two poly- 
nomials by splitting the division over addition and subtraction and by dividing 
out common powers of variables and greatest common divisors of the constant 
factors. Subsequently, fractions whose denominators have identical power prod- 
ucts are added. We used this optimization in the initial example by simplifying 
to ^ + The reader may check that without this simplification step we 
can only deduce that y < 43. 

To provide details of this optimization, given two monomials s and t, we 
denote by 



the result of performing this simplification operation on s and t. For example, 
[2^] equals ^, whereas [^] equals ^. 

In this approach we assume that the domains of the variables yi, . . .,yp-i, 
Up+i, . . ., t/n of mi do not contain 0. (One can easily show that this restriction is 
necessary here). For a monomial s involving variables ranging over the integer 
intervals that do not contain 0, the set int{s) either contains only positive num- 
bers or only negative numbers. In the first case we write sign{s) = + and in the 
second case we write sign{s) = — . 



The new domain of the variable Xj in the POLYNOMIAL INEQUALITY 
rule is defined using two sequences m'Q...m'^ and s'q...s'^ of extended arithmetic 
expressions such that 




Let S := {s[ \ i G [0..m] — {£}} and for an extended arithmetic expression t & S 
let It := {i e [0..m] — {£} \ s^ = t}. We denote then by pt the polynomial 
Sie/t "^^^ "^^^ domains are then defined by 




if sign{s) = +, and by 




if sign{s) = —. Here the int{s) notation used in the Correctness Lemma 1 is 
extended to expressions involving the division operator on real intervals in the 
obvious way. We define the int{.) operator applied to a bounded set of real 
numbers, as produced by the division and addition operators in the above two 
expressions for Dp to denote the smallest interval of real numbers containing 
that set. 

5 Constraint Propagation: Second Approach 

In this approach we limit our attention to a special type of polynomial con- 
straints, namely the ones of the form s op b, where s is a polynomial in which 
each variable occurs at most once and where b is an integer. We call such a con- 
straint a simple polynomial constraint. By introducing auxiliary variables that 
are equated with appropriate monomials we can rewrite each polynomial con- 
straint into a sequence of simple polynomial constraints. This allows us also to 
compute the integer interval domains of the auxiliary variable from the integer 
interval domains of the original variables. We apply then to the simple polyno- 
mial constraints the rules introduced in the previous section. 

To see that the restriction to simple polynomial constraints can make a dif- 
ference consider the constraint 

IQQx -y-lOy- z = 212 

in presence of the ranges x,y,z G [1..9]. We rewrite it into the sequence 

u = X ■ y, V = y ■ z, lOOu — lOv = 212 



where u,v are auxiliary variables, each with the domain [1..81]. 



It is easy to check that the POLYNOMIAL EQUALITY rule introduced 
in the previous section docs not yield any domain reduction when applied to 
the original constraint lOOx ■ y — lOy ■ z = 212. In presence of the discussed 
optimization the domain of x gets reduced to [1..3]. 

However, if wc repeatedly apply the POLYNOMIAL EQUALITY rule to 
the simple polynomial constraint lOOu — lOv = 212, we eventually reduce the 
domain of u to the empty set (since this constraint has no integer solution in the 
ranges u,v G [1--81]) and consequently can conclude that the original constraint 
lOOx ■ y — lOy ■ z — 212 has no solution in the ranges x,y,z € [1..9], without 
performing any search. 

6 Constraint Propagation: Third Approach 

In this approach wc focus on a small set of 'atomic' arithmetic constraints. We 
call an arithmetic constraint atomic if it is in one of the following two forms: 

— a linear constraint, 

— X ■ y = z. 

It is easy to see that using appropriate transformation rules involving aux- 
iliary variables we can transform each arithmetic constraint to a sequence of 
atomic arithmetic constraints. In this transformation, as in the second approach, 
the auxiliary variables are equated with monomials so we can easily compute 
their domains. 

The transformation to atomic constraints can strengthen the reduction. Con- 
sider for example the constraint 

u-x-y + l = v- x- y 

and ranges u e [1--2], v € [3.. 4], and x,y € [1--4]. The first approach without 
optimization and the second approach cannot find a solution without search. 
If, as a first step in transforming this constraint into a linear constraint, we 
introduce an auxiliary variable w to replace x ■ y, we are effectively solving the 
constraint 

U ■ W + 1 = V ■ w 

with the additional range w G [1..16], resulting in only one duplicate occurrence 
of a variable instead of two. With variable w introduced (or using the optimized 
version of the first approach) constraint propagation alone finds the solution 
u = 2, w = 3, X = 1, y = 1. 

We explained already in Section 3 how to reason about linear constraints. 
(We omitted there the treatment of the disequalities which is routine.) Next, we 
focus on the reasoning for the multiplication constraint x ■ y = z in presence of 
the non-empty ranges x € D^, y € Dy and z G D^. To this end we introduce the 
following three domain reduction rules: 



MULTIPLICATION 1 

{x-y = z\ X e Dx,y e Dy,z e D:,) 

{x-y = z; X e Dx, y e Dy, z e D^^n int{Dx ■ Dy)) 



MULTIPLICATION 2 

{x-y = z\ X € Dx,y& Dy,z e Dz) 

{x -y = z ; X G D^n int{Dz/Dy),y e Dy, z e D^) 



MULTIPLICATION 3 

{x-y = z; X e Dx,y e Dy,z e Dz) 

{x-y = z ; X e D^, yeDyD int{Dz/Dx), z e D^) 

The way we defined the multiphcation and the division of the integer intervals 
ensures that the MULTIPLICATION rules 1,2 and 3 are equivalence preserving. 
Consider for example the MULTIPLICATION 2 rule. Take some aeD^^heDy 
and c£ Dz such that a-h = c. Then a € {x & Z \ 3z ^ Dz^y £ Dy x ■ y = z}, so 
a e Dz/Dy and a fortiori a G int{Dz/Dy). Consequently a G D^d int{Dz/ Dy). 
This shows that the MULTIPLICATION 2 rule is equivalence preserving. 

The following example shows an interaction between all three MULTIPLI- 
CATION rules. 

Example 1. Consider the CSP 

{x-y = z; X e [1..20], y € [9..11], .z e [155.. 161]). (4) 

To facilitate the reading we underline the modified domains. An application 
of the MULTIPLICATION 2 rule yields 

{x-y = z; x& [16..16] , y e [9.. 11], z e [155.. 161]) 

since, as already noted in in Subsection 2.2, [155..161]/[9..11]) = [16. .16], and 
[1..20] n int([16..16]) = [16.. 16]. Applying now the MULTIPLICATION 3 rule 
we obtain 

{x-y = z; X G [16.. 16], 2/ G [10.. 10] , z G [155. .161]) 

since [155..161]/[16..16] = [10.. 10] and [9..11] n mt([10..10]) = [10.. 10]. Next, by 
the appUcation of the MULTIPLICATION 1 rule we obtain 

{x-y = z; X e [16..16],y G [10..10],z e [160.. 160] ) 

since [16. .16] • [10.. 10] = [160.. 160] and [155..161] n mt([160..160]) = [160.. 160]. 
So using all three multiplication rules we could solve the CSP (4). □ 



Now let us clarify why we did not define the division of the sets of integers 
Z and Y by 

ZjY := {z/y e Z \ y e Y, z e Z,y ^ 0}. 

The reason is that in that case for any set of integers Z we would have Z/ {0} = 0. 
Consequently, if we adopted this definition of the division of the integer inter- 
vals, the resulting MULTIPLICATION 2 and 3 rules would not be anymore 
equivalence preserving. Indeed, consider the CSP 

{x-y = z ; xe [-2..1], y e [0..0], ^ € [-8.. 10]). 

Then we would have [-8..10]/[0..0] = and consequently by the MULTIPLI- 
CATION 2 rule we could conclude 

{x-y = z; xGlD,ye [0..0], z e [-8.. 10]). 

So we reached an inconsistent CSP while the original CSP is consistent. 

In the remainder of the paper we will also consider variants of this third 
approach that allow squaring and exponentiation as atomic constraints. For this 
purpose we explain the reasoning for the constraint a; = in presence of the 
non-empty ranges x G and y G Dy, and for n > 1. To this end we introduce 
the following two rules in which to maintain the property that the domains are 
intervals we use the int{.) operation of Section 2: 

EXPONENTIA TION 
{x = y"' ] xeDx,ye Dy) 
{x = y^; xeD^n int{D^),y e Dy) 

ROOT EXTRACTION 

(x = ; X & Dx,y e Dy) 
{x ^y'"- ; X e D^,y e int{Dy n ^/DT)) 

To prove that these rules are equivalence preserving suppose that for some 
a G Dx and b G Dy wc have a = b". Then a S Dy, so a & int{Dy) and conse- 
quently a € DxC] int{Dy). Also 6 e \/Dx, so 6 e fl \/Dx, and consequently 
b e int{Dy n y/D^). 

7 Implementation Details 

In this section we describe the benchmark experiments that were performed 
to compare the proposed approaches. These experiments were performed using 
a single solver of the DICE (Distributed Constraint Environment) framework. 
DICE [10] is a framework for solver cooperation, implemented using techniques 
from coordination programming. It is developed around an experimental con- 
straint solver, called OpcnSolvcr, which is particularly suited for coordination. 
The coordination and cooperation aspects are irrelevant from the point of view 
of this paper. Relevant aspects of the OpenSolver are: 



— It implements a branch-and-infer tree search algorithm for constraint solv- 
ing. The inference stage corresponds to constraint propagation and is per- 
formed by repeated application of domain reduction functions (DRFs) that 
correspond to the domain reduction rules associated with the considered 
constraints. 

— This algorithm is abstract in the sense that the actual functionality is deter- 
mined by software plug-ins in a number of predefined categories. These cat- 
egories correspond to various aspects of the abstract branch-and-infer tree 
search algorithm. Relevant categories are: variable domain types, domain 
reduction functions, schedulers that control the application of the DRFs, 
branching strategies that split the search tree after constraint propagation 
has terminated, and several categories corresponding to different aspects of 
a search strategy that determine how to traverse a search tree. 

All experiments were performed using the Integerlnterval variable do- 
main type plug-in. Domains of this type consist of an indication of the type of 
the interval (bounded, unbounded, left/right-bounded, or empty), and a pair 
of arbitrary precision integer bounds. This plug-in, and the interval arithmetic 
operations on it are built using the GNU MP library [4]. 

The branching strategy that we used selects variables using the chronological 
ordering in which the auxiliary variables come last. The domain of the selected 
variable is split into two subdomains using bisection, so the resulting search trees 
are binary trees. In all experiments we searched for all solutions, traversing the 
entire search tree by means of depth-first leftmost-first chronological backtrack- 
ing. 

For the experiments in this paper a DRF plug-in has been developed that 
implements the domain reduction rules discussed in the previous sections. The 
scheduler plug-in used in the benchmarks keeps cycling through the sequence of 
DRFs, applying DRFs that have been scheduled for execution. When a DRF is 
applied, and some variable domain is modified, all DRFs that depend on these 
changes are scheduled for execution, including possibly the one that has just been 
applied. The cycling stops when no more DRFs are scheduled for execution, or 
when the domain of a variable becomes empty. 

As an alternative to cycling, the scheduler can be supplied with a schedule: 
a sequence of indices into the sequence of DRFs. The scheduler will then cycle 
through this schedule instead, and consider DRFs for application in the specified 
order. This is used in combination with the second and third approach, where we 
distinguish user constraints from the constraints that are introduced to define 
the values of auxiliary variables. Before considering for execution a DRF / that 
is part of the implementation of a user constraint, we make sure that all auxiliary 
variables that / relies on are updated. For this purpose, the indices of the DRFs 
that update these variables precede the index of / in the schedule. If / can 
change the value of an auxiliary variable, its index is followed by the indices of 
the DRFs that propagate back these changes to the variables that define the 
value of this auxiliary variable. 



For the third approach, there can be hierarchical dependencies between aux- 
ihary variables. Much like the HC4 algorithm of [2], the schedule specifies a 
bottom-up traversal of this hierarchy in a forward evaluation phase and a top- 
down traversal in a backward propagation phase before and after applying a 
DRF of a user constraint, respectively. In the forward evaluation phase, the 
DRFs that are executed correspond to rules MULTIPLICATION 1 and EXPO- 
NENTIATION. The DRFs of the backward propagation phase correspond to 
MULTIPLICATION 2 and 3, and ROOT EXTRACTION. It is easy to con- 
struct examples showing that the use of hierarchical schedules can be beneficial 
compared to cycling through the rules. 

The proposed approaches were implemented by first rewriting arithmetic 
constraints to polynomial constraints, and then to a sequence of DRFs that 
correspond with the rules of the approach used. We considered the following 
methods: 

la the first approach, discussed in Section 4, 

lb the optimization of the first approach discussed at the end of Section 4 that 
involves dividing out common powers of variables, 

2a the second approach, discussed in Section 5. The conversion to simple poly- 
nomial constraints is implemented by introducing an auxiliary variable for 
every non-linear monomial. This procedure may introduce more auxiliary 
variables than necessary. 

2b an optimized version of approach 2a, where we stop introducing auxiliary 
variables as soon as the constraints contain no more duplicate occurrences 
of variables. 

3a the third approach, discussed in Section 6, allowing only linear constraints 

and multiplication as atomic constraints. 
3b idem, but also allowing x = as an atomic constraint. 
3c idem, allowing x = for all n > 1 as an atomic constraint. 

Approaches 2 and 3 involve an extra rewrite step, where the auxiliary vari- 
ables are introduced. The resulting CSP is then rewritten according to approach 
la. During the first rewrite step the hierarchical relations between the auxiliary 
variables are recorded and the schedules are generated as a part of the second 
rewrite step. For approaches 2b and 3 the question of which auxiliary variables 
to introduce is an optimization problem in itself. Some choices result in more 
auxiliary variables than others. We have not treated this issue as an optimiza- 
tion problem but relied on heuristics. We arc confident that these yield a realistic 
implementation. In our experiments we used the following benchmarks. 

Cubes The problem is to find all natural numbers n < 1000 that are a sum of 
four different cubes, for example 

1^ + 2^ + 'd'-'' + 4'' = 100. 

This problem is modeled as follows: 

(1 < Xi, Xi < X2 — 1, X2 < X3 — 1, X3 < X4 — 1, X4 < 71, 

Xi + X2 + x^ + X4 = n; n G [1..1000], xi,X2,xs,Xi e Z) 



Opt We are interested in finding a solution to the constraint + = in the 
integer interval [1.. 100000] for which the value o{2x ■ y — z is maximal. 



Fractions This problem is taken from [9]: find distinct nonzero digits such that 
the following equation holds: 

A D G _ 

There is a variable for each letter. The initial domains are [1..9]. To avoid sym- 
metric solutions an ordering is imposed: 

A D G 

'bc - ^ - m 

Also two redundant constraints are added: 

3A>i and 3|^<1 

Because division is not present in the arithmetic expressions, the above con- 
straints are multiplied by the denominators of the fractions to obtain arithmetic 

constraints. 

Two representations for this problem were studied: 

— fractionsl in which five constraints arc used: one equality and four inequal- 
ities for the ordering and the redundant constraints, 

— fractions2, used in [9] , in which three auxiliary variables, BG, EF and HI, 
are introduced to simplify the arithmetic constraints: EC = lOB + G, EF = 
lOE + F, and HI = lOH + 1. 

Additionally, in both representations, 36 disequalities A ^ B, A ^ G, H ^ I 
are used. 



Kyoto The problem is to find the number n such that the alphanumeric equation 

KYOTO 
KYOTO 
+KYOTO 
TOKYO 

has a solution in the base-n number system. Our model uses a variable for each 
letter and one variable for the base number. The variables K and T may not 
be zero. There is one large constraint for the addition, 6 disequalities K ^ Y 
... T ^ O and four constraints stating that the individual digits K,Y,0,T, 
arc smaller than the base number. To spend some CPU time, we searched base 
numbers 2. .100. 

^ V. Dubrovsky and A. Shvetsov. Quantum cyberteaser: May/June 1995, 
http: / / www.nsta.org/ quantum/kyotoarc.asp 



8 Results 



Tables 1 and 2 compare the proposed approaches on the problems defined in the 
previous section. The first two columns of table 1 list the number of variables and 
the DRFs that were used. Column nodes lists the size of the search tree, including 
failures and solutions. The next two columns list the number of times that a 
DRF was executed, and the percentage of these activations that the domain of a 
variable was actually modified. For the opt problem, the DRF that implements 
the optimization is not counted, and its activation is not taken into account. The 
elapsed times in the last column are the minimum times (in seconds) recorded 
for 5 runs on a 1200 MHz Athlon CPU. 
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16 


59 
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3,746,532 


23.26 


10.466 







Table 1. Statistics and comparison with other solvers 



Table 2 lists measured numbers of basic interval operations. Note that for 
approach lb, there are two versions of the division and addition operations: 
one for integer intervals, and one for intervals of reals of which the bounds are 
rational numbers (marked Q). Columns multl and multF list the numbers of 



multiplications of two integer intervals, and of an integer interval and an integer 
factor, respectively. These are different operations in our implementation. 
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1 


1 
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5,324 


7,504 


14,868 



Table 2. Measured numbers (thousands) of interval operations 



For the cubes and opt problems, the constraints are already in simple form, so 
approaches la, lb and 2b are identical. Also all non-linear terms involve either a 
multiplication or an exponentiation, so also approaches 2a and 3c arc the same. 
The results of these experiments clearly show the disadvantage of implementing 
exponentiation by means of multiplication: the search space grows because we 
increase the number of variable occurrences and lose the information that it is 
the same number that is being multiplied. For opt and approach 3a, the run did 
not complete within reasonable time and was aborted. 



Columns E and I of table 1 compare the propagation achieved by our ap- 
proaches with two other systems, respectively ECL^PS*^ Version 5.6^ using the 
ic library, and ILOG Solver 5.1^ using type ILOINT. For this purpose we ran 
the test problems without search, and compared the results of constraint prop- 
agation. A mark '=' means that the computed domains arc the same, '+' that 
our approach achieved stronger propagation than the solver that we compare 
with, and '-' that propagation is weaker. For cubes, ECL'PS* computes the same 
domains as those computed according to approach 3b, so here the reduction is 
stronger than for 3a, but weaker than for the other approaches. For opt ECL'PS'^ 
and ILOG Solver compute the same domains. These domains are narrower than 
those computed according to approaches 3a and 3b, but the other approaches 
achieve stronger reduction. In all other cases except for kyoto and approach lb 
the results of all three solvers are the same. 

For both representations for the fractions puzzle, the symbolic manipulation 
of approach lb is able to achieve a significant reduction of the search tree, but 
this is not reflected in the timings. For fractionsl the elapsed time even increases. 
The reason is that computing the domain updates involves adding intervals of 
real numbers. The arithmetic operations on such intervals are more expensive 
than their counterparts on integer intervals, because the bounds have to be main- 
tained as rational numbers. Arithmetic operations on rational numbers are more 
expensive because they involve the computation of greatest common divisors. 
For kyoto the symbolic manipulation did not reduce the size of the search tree, 
so the effect is even more severe. 

In general, the introduction of auxiliary variables leads to a reduction of 
the number of interval operations compared to approach la. The reason is that 
auxiliary variables prevent the evaluation of subexpressions that did not change. 
This effect is strongest for fractionsl, where the main constraint contains a 
large number of different power products. Without auxiliary variables all power 
products are evaluated for every POLYNOMIAL EQUALITY rule defined by 
this constraint, even those power products the variable domains of which did 
not change. With auxiliary variables the intervals for such unmodified terms are 
available immediately, which leads to a significant reduction of the number of 
interval multiplications. 

The efi^ect that stronger reduction is achieved as a result of introducing aux- 
iliary variables, mentioned in Section 6, is seen for both representations of the 
fractions benchmark. The effect described in Section 5 is not demonstrated by 
these experiments. 

If we don't consider the symbolic manipulation of approach lb, approach 
3c leads to the smallest total number of interval operations in all cases, but 
the scheduling mechanism discussed in Section 7 is essential for a consistent 
good performance. If for example the schedule is omitted for opt, the number 



^ ECL'PS'^ Constraint Logic Programming System. See http://www- 

icparc.doc.ic.ac.uk/eclipse 
^ See http://www.ilog.com 



of interval operations almost triples, and performance of approach 2a and 3c is 
then much worse than that of approach la. 

The total numbers of interval operations in table 2 do not fully explain all 
differences in elapsed times. One of the reasons is that different interval opera- 
tions have different costs. Especially the preprocessing of the numerator interval 
for integer interval division, discussed in Subsection 2.2, is potentially expen- 
sive, which may explain why for opt, approach la runs faster than approach 2a, 
even though the total number of interval operations is higher. Among the many 
other factors that may be of influence, some overhead is involved in applying a 
DRF, so if the number of applications differs significantly for two experiments, 
this probably influences the elapsed times as well {cubes, la, 2a, opt, la, 2a, 
fractions2, 2a, 2b). The elapsed times are not the only measure that is subject 
to implementation details. For example, we implemented division by a constant 
interval [—1.. — 1] as multiplication by a constant, which is more efficient in our 
implementation. Such decisions are reflected in the numbers reported in table 2. 

9 Discussion 

In this paper we discussed a number of approaches to constraint propagation 
for arithmetic constraints on integer intervals. To assess them we implemented 
them using the DICE (Distributed Constraint Environment) framework of [10], 
and compared their performance on a number of benchmark problems. We can 
conclude that: 

— Implementation of exponentiation by multiplication gives weak reduction. 
In our third approach x = y" should be an atomic constraint. 

— The optimization of the first approach, where common powers of variables 
arc divided out, can significantly reduce the size of the search tree, but 
the resulting reduction steps rely heavily on the division and addition of 
rational numbers. These operations can be expected to be more expensive 
than their integer counterparts, because they involve the computation of 
greatest common divisors. 

— Introducing auxiliary variables can be beneficial in two ways: it may strength- 
en the propagation, as discussed in Sections 5 and 6, and it may prevent the 
evaluation of subexpressions the variable domains of which did not change. 

— As a result, given a proper scheduling of the rules, the second and third 
approach perform better than the first approach without the optimization, 
in terms of numbers of interval operations. Actual performance depends on 
many implementation aspects. However for our test problems the results of 
variants 2a, 2b and 3c do not diff'er significantly. 

In general, our implementation is slow compared to, for example, ILOG 

Solver. A likely cause is that we use arbitrary precision integers. Wc chose this 
representation to avoid having to deal with overfiow, but an additional benefit 
is that large numbers can be represented exactly. 



A different approach would be to use floating-point arithmetic and then round 
intervals inwards to the largest enclosed integer interval. This was suggested in [3] 
and implemented in for example RealPaver^. A benefit of this inward rounding 
approach is that all algorithms that were developed for constraints on the reals 
are immediately available. A disadvantage is that for large numbers no precise 
representation exists, i.e., the interval defined by two consecutive floating-point 
numbers contains more than one integer. But it is debatable whether an exact 
representation is required for such large numbers. 

We realize that the current set of test problems is rather limited. In addition 
to puzzles, some more complex non-linear integer optimization problems should 
be studied. We plan to further evaluate the proposed approaches on non-linear 
integer models for the SAT problem. Also we would like to study the relationship 
with the local consistency notions that have been defined for constraints on the 
reals and give a proper characterization of the local consistencies computed by 
our reduction rules. 

Note This work was performed during the first author's stay at the School of 
Computing of the National University of Singapore. The work of the second 
author was supported by NWO, The Netherlands Organization for Scientific 
Research, under project number 612.069.003. 
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